REPORT  DOCUMENTATION  PAGE 


Form  Approved  OMB  NO.  0704-0188 


The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments 
regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggesstions  for  reducing  this  burden,  to  Washington 
Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington  VA,  22202-4302. 
Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  oenalty  for  failing  to  comply  with  a  collection 
of  information  if  it  does  not  display  a  currently  valid  OMB  control  number. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 


2.  REPORT  TYPE 

Final  Report 


1.  REPORT  DATE  (DD-MM-YYYY) 

07-01-2016 


4.  TITLE  AND  SUBTITLE  5a.  CONTRACT  NUMBER 

Final  Report:  Fast,  Multiscale  Algorithms  for  Wave  Propagation  W91  INF-09- 1-0344 
in  Heterogeneous  Environments  5b.  GRANT  NUMBER 


3.  DATES  COVERED  (From  -  To) 

16-Jun-2009-  15-Jan-2015 


6.  AUTHORS 
Thomas  Hagstrom 


5c.  PROGRAM  ELEMENT  NUMBER 
611102 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


7.  PERFORMING  ORGANIZATION  NAMES  AND  ADDRESSES 

Southern  Methodist  University 

Research  Administration 

6425  Boaz  Lane,  Room  G05 

Dallas,  TX  75275  -0302 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS 
(ES) 

U.S.  Army  Research  Office 
P.O.Box  12211 

Research  Triangle  Park,  NC  27709-2211 


8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 
ARO 


11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

56543-MA.38 


12.  DISTRIBUTION  AVAIL1B1LITY  STATEMENT 
Approved  for  Public  Release;  Distribution  Unlimited 


13.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  should  not  contrued  as  an  official  Department 
of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


14.  ABSTRACT 

The  objective  of  this  research  project  was  to  further  develop  and  integrate  numerical  methods  for  the  fast  and 
accurate  simulation  of  wave  propagation  problems  in  the  time  domain.  In  support  of  the  long-term  goal  of  creating 
high-quality  software  for  simulating  waves,  we  seek  methods  which  are  not  only  efficient,  but  which  are  reliable  in 
that  both  their  stability  and  the  accuracy  of  the  results  are  essentially  guaranteed.  In  support  of  this  goal  we  have 
developed:  (i.)  convenient  implementations  of  optimal  local  radiation  boundary  sequences  for  isotropic  waves,  with 

-a  m  1  n  i  vt  ,  t  4  t  t  r»  /%4r  r  /%  T  m  i!  4  a  /-a  r—r  s~a  n  n  Tr\  «•  A  Tri-^rTT  rrtl  1  In  n  •  /i  «  \  /-virt-  n  ■*  n  /-v  T 


15.  SUBJECT  TERMS 
numerical  methods,  wave  propagation 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

15.  NUMBER 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

ABSTRACT 

OF  PAGES 

UU 

UU 

UU 

UU 

Thomas  Hagstrom 


19b.  TELEPHONE  NUMBER 
214-768-4338 


Standard  Form  298  (Rev  8/98) 
Prescribed  by  ANSI  Std.  Z39.18 


Report  Title 

Final  Report:  Fast,  Multiscale  Algorithms  for  Wave  Propagation  in  Heterogeneous  Environments 

ABSTRACT 

The  objective  of  this  research  project  was  to  further  develop  and  integrate  numerical  methods  for  the  fast  and  accurate  simulation  of  wave 
propagation  problems  in  the  time  domain.  In  support  of  the  long-term  goal  of  creating  high-quality  software  for  simulating  waves,  we  seek 
methods  which  are  not  only  efficient,  but  which  are  reliable  in  that  both  their  stability  and  the  accuracy  of  the  results  are  essentially 
guaranteed.  In  support  of  this  goal  we  have  developed:  (i.)  convenient  implementations  of  optimal  local  radiation  boundary  sequences  for 
isotropic  waves,  with  implementations  in  a  wide  variety  of  popular  discretization  schemes  for  Maxwell's  equations;  (ii.)  extensions  of  these 
sequences  to  more  complex  systems  arising  in  linear  elasticity;  (iii.)  new  highly  efficient  energy-stable  discretization  schemes  on  structured 
grids  -  these  include  methods  based  on  Hermite  interpolation  and  compact  difference  schemes  constructed  using  Galerkin  techniques;  (iv.) 
stable  coupling  of  the  efficient  structured  grid  methods  with  upwind  discontinuous  Galerkin  methods  defined  on  unstructured  grids  -  using 
hybrid  grids  allows  us  to  treat  very  complex  geometry  with  efficency  comparable  to  simple  domains;  (v.)  natural  upwind  discontinuous 
Galerkin  discretizations  for  wave  equations  in  second  order  form  -  using  the  second  order  form  for  complex  systems  results  in  fewer 
dependent  variables. 


Enter  List  of  papers  submitted  or  published  that  acknowledge  ARO  support  from  the  start  of 
the  project  to  the  date  of  this  printing.  List  the  papers,  including  journal  references,  in  the 
following  categories: 

(a)  Papers  published  in  peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


01/04/2016  21.00 

01/04/2016  23.00 

01/04/2016  24.00 

01/04/2016  25.00 

01/04/2016  27.00 

08/31/2012  1.00 

08/31/2012  11.00 

08/31/2012  10.00 

08/31/2012  7.00 

08/31/2012  2.00 

08/31/2012  3.00 

08/31/2012  4.00 

09/17/2013  16.00 


M.  Beroiz,  T.  Hagstrom,  S.R.  Lau,  R.H.  Price.  Multidomain,  sparse,  spectral-tau  method  for  helically 
symmetric  flow, 

Computers  &  Fluids,  (10  2014):  0.  doi:  10.1016/j.compfluid.2014.05.028 

Xi  (Ronald)  Chen,  Daniel  Appelo,  Thomas  Hagstrom.  A  hybrid  Hermite-discontinuous  Galerkin  method 
for  hyperbolic  systems  with  application  to  Maxwell?s  equations, 

Journal  of  Computational  Physics,  (01  2014):  501.  doi:  10.101 6/j.jcp.201 3.09.046 

Thomas  Hagstrom.  High-resolution  difference  methods  with  exact  evolution  for  multidimensional  waves, 
Applied  Numerical  Mathematics,  (07  2015):  114.  doi:  10.1016/j.apnum.2014.07.001 

Daniel  Baffet,  Thomas  Hagstrom,  Dan  Givoli.  Double  Absorbing  Boundary  Formulations  for  Acoustics  and 
Elastodynamics, 

SIAM  Journal  on  Scientific  Computing,  (01  2014):  1277.  doi:  10.1137/130928728 

Daniel  Appelo,  Thomas  Hagstrom.  A  new  discontinuous  Galerkin  formulation  for  wave  equations  in 
second  order  form, 

SIAM  J  Numer  Anal,  (12  2015):  2705.  doi: 

Ronald  Chen,  Thomas  Hagstrom.  P-adaptive  Hermite  methods  for  initial  value  problems, 

Mathematical  Modeling  and  Numerical  Analysis,  (05  2012):  545.  doi: 

Thomas  Hagstrom.  HIGH-ORDER  RADIATION  BOUNDARY  CONDITIONS  FORSTRATIFIED  MEDIA 
AND  CURVILINEAR  COORDINATES, 

Journal  of  Computational  Acoustics,  (06  2012):  0.  doi: 

Daniel  Baffet,  Jacobo  Bielak,  Dan  Givoli,  Thomas  Hagstrom,  Daniel  Rabinovich.  Long-time  stable  high- 
order  absorbing  boundary  conditions  for  elastodynamics, 

Computer  Methods  in  Applied  Mechanics  and  Engineering,  (10  2012):  20.  doi: 

Daniel  Appelo,  Thomas  Hagstrom.  On  Advection  by  Hermite  Methods, 

Pacific  Journal  of  Applied  Mathematics,  (05  2012):  0.  doi: 

Thomas  Hagstrom,  George  Hagstrom.  Grid-stabilization  of  high-order  one-sided  differencing  II.  Second- 
orderwave  equations, 

Journal  of  Computational  Physics,  (01  2012):  0.  doi: 

Eliane  Becache,  Dan  Givoli,  Thomas  Hagstrom,  Kurt  Stein.  Complete  radiation  boundary  conditions  for 
convective  waves, 

COMMUNICATIONS  IN  COMPUTATIONAL  PHYSICS,  (10  2011):  610.  doi: 

Daniel  Rabinovich,  Dan  Givoli,  Jacobo  Bielak,  Thomas  Hagstrom.  A  finite-element  scheme  with  a  high- 
order  absorbing  boundary  condition  for  elastodynamics, 

Computer  Methods  in  Applied  Mechanics  and  Engineering,  (06  2011):  2048.  doi: 

Daniel  Baffet,  Jacobo  Bielak,  Dan  Givoli,  Thomas  Hagstrom,  Daniel  Rabinovich.  Long-time  stable  high- 
order  absorbing  boundary  conditions  for  elastodynamics, 

Computer  Methods  in  Applied  Mechanics  and  Engineering,  (05  2012):  241.  doi: 


09/17/2013  18.00  Daniel  Rabinovich,  Dan  Givoli,  Thomas  Hagstrom,  Jacobo  Bielak  .  STRESS-VELOCITY  COMPLETE 
RADIATIONBOUNDARY  CONDITIONS, 

Journal  of  Computational  Acoustics,  (04  2013):  1350003.  doi: 

TOTAL:  14 


Number  of  Papers  published  in  peer-reviewed  journals: 


(b)  Papers  published  in  non-peer-reviewed  journals  (N/A  for  none) 


Received  Paper 


01/05/2016  31.00  Leslie  Greengard,  Thomas  Hagstrom,  Shidong  Jiang  .  Extension  of  the  Lorenz-Mie-Debye  method 
forelectromagnetic  scattering  to  the  time-domain, 

Journal  of  Computational  Physics,  (07  2015):  98.  doi: 

01/05/2016  33.00  Daniel  Rabinovich,  Dan  Givoli,  Jacobo  Bielak,  Thomas  Hagstrom.  The  double  absorbing  boundary 
method  forel  astodynamics  in  homogeneous  and  layered  media, 

Advanced  Modeling  and  Simulation  in  Engineering  Sciences,  (04  2015):  0.  doi: 

TOTAL:  2 


Number  of  Papers  published  in  non  peer-reviewed  journals: 


(c)  Presentations 


"A  few  things  I  learned  from  Herb  Keller",  Keller/Isaacson  Memorial  Conference  on  Numerical  Analysis,  Courant  Institute,  New  York, 
2009. 

"An  Introduction  to  Radiation  Boundary  Conditions  for  Time-Domain  Scattering  Problems",  IMA  PI  Summer  School  on  Computational 
Wave  Propagation,  Michigan  State  University,  2010. 

"Complete  Radiation  Boundary  Conditions:  Theory  and  Application",  SIAM  Annual  Meeting,  Pittsburgh,  2010. 

"Complete  radiation  boundary  conditions  on  rectangular  domains",  AMS  Midwest  Sectional  Meeting,  Notre  Dame,  2010. 

"Approximate  Radiation  Boundary  Conditions  for  Time-Dependent  Waves",  Applied  Inverse  Problems,  Texas  A&M,  2011. 

"Towards  the  Ultimate  Solver  for  the  Wave  Equation  in  the  Time  Domain",  Applied  Math  and  Imaging  Sciences  Workshop,  UT  Pan- 
American,  2011. 

"Complete  radiation  boundary  conditions",  Cambridge  Workshop  on  Computational  Wave  Propagation,  Isaac  Newton  Institute, 
Cambridge,  UK  2010. 

"HP -Refinement  and  DG  Coupling  for  Hermite  Methods"  2011  Finite  Element  Rodeo,  Texas  A&M. 

"Complete  Wave  Representations  and  Optimal  Local  Radiation  Boundary  Conditions",  SIAM  Conference  on  Comptational  Science  and 
Engineering,  Reno,  2011. 

"Accurate  Methods  for  Time-Domain  Scattering",  SIAM  Conference  on  Comptational  Science  and  Engineering,  Reno,  2011. 

"Boundary  conditions  for  simulating  waves",  Eigenvalues/singular  values  and  fast  PDE  algorithms:  acceleration,  conditioning,  and 
stability,  Banff  International  Research  Station,  Banff,  CA  2012. 

"Discretization  methods  for  waves".  Nonlinear  solvers  for  high-intensity  focused  ultrasound  with  application  to  cancer  treatment, 

AIMS,  Palo  Alto,  2012. 

"Hermite  methods  for  aeroacoustics",  Acoustics  2012,  Hong  Kong. 

"Hermite  spectral  elements",  Finite  Element  Rodeo  2012,  Houston. 

"Hermite  methods  for  hyperbolic  systems:  Basic  theory",  SIAM  Conference  on  Computational  Science  and  Engineering,  Boston  2013 
"Hybrid  Hermite-DG  Methods  for  Hyperbolic  IB  VPs",  Finite  Element  Rodeo,  LSU  2013 

"Towards  the  Ultimate  Solver  for  Wave  Equations  in  the  Time  Domain",  Forum  on  Scientific  and  Engineering  Computing,  Chinese 
Academy  of  Sciences,  Beijing  2013 

"Boundary  conditions  for  high-resolution  simulations  of  waves",  International  Conference  on  Difference  Schemes  and  their  Applications, 
Moscow  2013 

"Adaptive  and  hybridized  Hermite  methods  for  initial-boundary  value  problems",  MAFELAP  2013,  London 

"High-Order  Energy-Stable  Methods  for  Hyperbolic  IB  VPs  on  Structured-Unstructured  Grids",  Finite  Element  Rodeo,  Austin,  TX  2014 
"Dispersion  and  Dissipation  of  Hermite  Methods  in  Multiple  Dimensions",  C.-Y.  Jang,  Finite  Element  Rodeo,  Austin,  TX  2014 
"High-resolution  upwind  methods  for  waves",  ICOSAHOM,  Salt  Lake  City  2014 
"Robust  high-order  methods  for  waves",  ICOSAHOM,  Salt  Lake  City  2014  (Plenary  Talk) 

"Energy  stable  difference  approximations  to  double  absorbing  boundary  conditions",  K.  Juhnke,  ICOSAHOM,  Salt  Lake  City  2014 

"Multidimensional  dissipation  and  dispersion  analysis  and  Lp  stability  for  Hermite  methods",  C.-Y.  Jang,  ICOSAHOM,  Salt  Lake  City 
2014 


"Hybrid  Methods  for  Elastic  Waves",  World  Congress  on  Computational  Mechanics,  Barcelona  2014 


Number  of  Presentations:  26.00 


Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


Received  Paper 


TOTAL: 


Number  of  Non  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 

Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


Received  Paper 


01/05/2016  29.00  Thomas  Hagstrom,  Jeffrey  Banks  .  Galerkin  Difference  Methods  from  Bandlimited  Interpolation  Functions, 
Waves  2015.  20-JUL-15,  .  :  , 

01/05/2016  37.00  Thomas  Hagstrom.  Radiation  boundary  conditions  for  time-domain  scattering  problems, 

Oberwohlfach  Workshop  on  Computational  Electromagnetism  and  Acoustics.  14-FEB-10,  .  :  , 

01/05/2016  36.00  Fritz  Juhnke,  Thomas  Hagstrom.  A  DG  discretization  of  a  Double  Absorbing  Boundary, 

Waves  2015.  20-JUL-1 5,  .:  , 

01/05/2016  35.00  John  Lagrone,  Thomas  Hagstrom.  Double  Absorbing  Boundaries  for  Finite  Difference  Time  Domain 
Electromagnetics, 

Waves  2015.  20-JUL-15,  .  :  , 

08/31/2012  5.00  Daniel  Appelo,  Matthew  Inkman,  Thomas  Hagstrom,  Tim  Colonius.  Recent  progress  on  Hermite  methods 
for  aeroacoustics, 

17th  AIAA/CEAS  Aeroacoustics  Conference  .  08-JUN-11,  .  :  , 

08/31/2012  6.00  Thomas  Hagstrom,  Daniel  Appelo,  Chang  Young  Jang  .  HERMITE  METHODS  FOR  HYPERBOLIC- 
PARABOLIC  SYSTEMS, 

Waves  201 1. 25-JUL-11,  .  :  , 

08/31/2012  8.00  Chang  Young  Jang,  Daniel  Appelo,  Tim  Colonius,  Thomas  Hagstrom,  Matthew  Inkman.  An  Analysis  of 
Dispersion  and  Dissipation  Properties  of  Hermite  Methods  and  its  Application  to  Direct  Numerical 
Simulation  of  Jet  Noise, 

18th  AIAA/CEAS  Aeroacoustics  Conference.  04-JUN-12,  .  :  , 

08/31/2012  12.00  Kurt  Stein,  Thomas  Hagstrom.  Complete  Radiation  Boundary  Conditions:  Corners  and  Edges, 

Waves  201 1. 25-JUL-11,  .:  , 


TOTAL: 


8 


Number  of  Peer-Reviewed  Conference  Proceeding  publications  (other  than  abstracts): 


(d)  Manuscripts 


Received  Paper 


01/05/2016  30.00  Daniel  Appelo,  Thomas  Hagstrom.  An  energy-based  discontinuous  Galerkin  discretizationof  the  elastic 
wave  equation  in  second  order  form, 

Comput.  Meth.  Appl.  Mech.  Eng  (07  2015) 

01/05/2016  34.00  Charles  Epstein,  Leslie  Greengard,  Thomas  Hagstrom  .  On  the  stability  of  time-domain  integral  equationsf 
or  acoustic  wave  propagation, 

Discrete  and  Continuous  Dynamical  Systems  (08  2015) 

01/05/2016  32.00  Jeffrey  Banks,  Thomas  Hagstrom  .  On  Galerkin  Finite  Difference  Methods, 

Journal  of  Computational  Physics  (08  2015) 

08/31/2012  13.00  Seungil  Kim,  Thomas  Hagstrom.  COMPLETE  RADIATION  BOUNDARY  CONDITIONS  FOR 
THEHELMHOLTZ  EQUATION  I:  WAVEGUIDES, 

Mathematics  of  Computation  (06  2012) 

09/17/2013  14.00  Daniel  Baffet,  Thomas  Hagstrom,  Dan  Givoli.  Double  Absorbing  Boundary  Formulations  for  Acoustics  and 
Elastodynamics, 

SIAM  J  Scientific  Computing  (08  2013) 

09/17/2013  15.00  Xi  Chen,  Daniel  Appelo,  Thomas  Hagstrom.  A  Hybrid  Hermite  -  Discontinuous  Galerkin  Method 
forHyperbolic  Systems  with  Application  to  Maxwell’s  Equations, 

Journal  of  Computational  Physics  (06  2013) 

09/17/2013  20.00  Leslie  Greengard,  Thomas  Hagstrom,  Shidong  Jiang.  The  solution  of  the  scalar  wave  equation  in  the 
exterior  of  a  sphere, 

Journal  of  Computational  Physics  (08  2013) 

09/17/2013  19.00  Thomas  Hagstrom,  Dan  Givoli,  Daniel  Rabinovich,  Jacobo  Bielak.  The  Double  Absorbing  Boundary 
Method, 

Journal  of  Computational  Physics  (05  2013) 

TOTAL:  8 


Number  of  Manuscripts: 

Books 


Received  Book 


TOTAL: 


Received 


Book  Chapter 


01/05/2016  28  00  Daniel  Appelo,  Thomas  Hagstrom  .  Solving  PDEs  with  Hermite  Interpolation,  Switzerland:  Springer-Verlag 
~  ,  (12  2015) 

TOTAL:  1 


Patents  Submitted 


Patents  Awarded 


Ford  Research  Fellow,  SMU 


Awards 


Graduate  Students 


NAME 

PERCENT  SUPPORTED  Discipline 

Kurt  Stein 

1.00 

Chang-Young  Jang 

0.12 

John  Lagrone 

0.25 

Fritz  Juhnke 

0.25 

FTE  Equivalent: 

1.62 

Total  Number: 

4 

Names  of  Post  Doctorates 


NAME 

PERCENT  SUPPORTED 

Seungil  Kim 

1.00 

FTE  Equivalent: 

1.00 

Total  Number: 

1 

Names  of  Faculty  Supported 


NAME 

PERCENT  SUPPORTED  National  Academy  Member 

Thomas  Hagstrom 

0.13 

FTE  Equivalent: 

0.13 

Total  Number: 

1 

Names  of  Under  Graduate  students  supported 

NAME 

FTE  Equivalent: 

Total  Number: 

PERCENT  SUPPORTED 

- 1 

Student  Metrics 

This  section  only  applies  to  graduating  undergraduates  supported  by  this  agreement  in  this  reporting  period 


The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period: .  0.00 

The  number  of  undergraduates  funded  by  this  agreement  who  graduated  during  this  period  with  a  degree  in 

science,  mathematics,  engineering,  or  technology  fields: . 0.00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  continue 

to  pursue  a  graduate  or  Ph.D.  degree  in  science,  mathematics,  engineering,  or  technology  fields: . 0.00 

Number  of  graduating  undergraduates  who  achieved  a  3.5  GPA  to  4.0  (4.0  max  scale): . 0.00 

Number  of  graduating  undergraduates  funded  by  a  DoD  funded  Center  of  Excellence  grant  for 

Education,  Research  and  Engineering: . o.OO 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  intend  to  work 

for  the  Department  of  Defense . 0.00 

The  number  of  undergraduates  funded  by  your  agreement  who  graduated  during  this  period  and  will  receive 

scholarships  or  fellowships  for  further  studies  in  science,  mathematics,  engineering  or  technology  fields: . 0.00 

Names  of  Personnel  receiving  masters  degrees 

NAME 
John  Lagrone 
Fritz  Juhnke 

Total  Number:  2 


Names  of  personnel  receiving  PHDs 

NAME 

Kurt  Stein 

Chang-Young  Jang 

Total  Number: 

2 

Names  of  other  research  staff 

NAME 

PERCENT  SUPPORTED 

FTE  Equivalent: 

Total  Number: 

Sub  Contractors  (DD882) 

Inventions  (DD882) 


Scientific  Progress 


Technology  Transfer 
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1  Introduction 

The  overall  goal  of  this  research  is  the  development  and  application  of  efficient  algorithms  for 
solving  wave  propagation  problems  in  the  time  domain,  with  a  particular  emphasis  on  problems 
in  electromagnetism,  elasticity,  and  aeroacoustics.  The  basic  challenges  in  building  such  algo¬ 
rithms  are  rooted  in  the  physical  nature  of  waves  themselves;  their  defining  feature  is  the  ability  to 
propagate  large  distances  relative  to  the  wavelength,  carrying  information  about  their  sources  and 
the  medium  through  which  they  have  traveled.  Thus  typical  wave  propagation  problems  exhibit 
multiple  spatiotemporal  scales. 

Two  crucial  components  of  the  highly-efficient,  general-purpose  wave  simulator  we  envision  are 

•  Reliable,  low-cost  methods  for  truncating  the  computational  domain  in  the  near  field,  thus 
avoiding  sampling  the  wave  field  throughout  space; 

•  Robust,  high-resolution  volume  discretizations  providing  guaranteed  accuracy  using  minimal 
degrees-of-freedom  per  wavelength. 

We  have  made  substantial  progress  on  both  of  these  issues,  as  we  will  detail  below. 

2  Radiation  Boundary  Conditions  and  Integral  Equations 

Our  recently-developed  theory  of  complete  radiation  boundary  conditions  (CRBCs)  solves  a  long¬ 
standing  problem  in  computational  wave  propagation;  namely  how  to  truncate  the  domain  in  the 
vicinity  of  regions  of  interest  in  such  a  way  that  reflections  from  the  computational  boundary  can 
be  efficiently  eliminated  to  any  desired  accuracy  [26,  22,  7,  23,  9].  CRBCs  are  the  provably 
optimal  method  for  achieving  this  in  the  case  of  isotropic  waves,  including  the  scalar  wave 
equation  and  Maxwell’s  equations.  In  particular,  as  proven  in  [26],  we  can  guarantee  an  accuracy,  e, 
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over  a  time  interval,  T,  assuming  a  separation,  5 ,  of  scatterers  and  sources  from  the  computational 
boundary  using  a  CRBC  with  p  auxiliary  variables  per  field  with 


p  oc  In 


•  In 


In  this  project  we  have  developed  full  implementations  of  CRBCs  for  both  exterior  and  waveg¬ 
uide  problems  in  three  space  dimensions,  devised  a  new  formulation  better  suited  to  second  order 
formulations  and  staggered- grid  codes,  extended  their  construction  to  elastic  waves  and  stratified 
media,  completed  an  open-source  software  library  to  enable  their  general  use.  In  addition  we  have 
initiated  a  study  of  numerical  methods  and  proper  formulations  of  time-domain  integral  equations 
for  scattering  problems. 


2.1  Implementation  of  CRBCs  with  Corners  and  Edges 


As  the  major  component  of  his  doctoral  thesis  in  computational  mathematics  at  SMU,  Kurt  Stein, 
supporated  as  an  RA  on  this  project,  completed  the  development  and  initial  testing  of  a  high- 
order,  parallelized,  3  +  1  dimensional  structured  grid  solver  for  first  order  hyperbolic  systems,  with 
an  implementation  of  CRBCs  built  in.  The  implementation  uses  grid-stabilized  finite  difference 
methods  [24,  25]  up  through  12th  order.  The  mathematical  structure  of  CRBCs  involves  the 
evolution  of  a  hyperbolic  system  of  auxiliary  variables  defined  only  on  the  domain  boundary.  For 
example,  along  a  boundary  with  normal  in  the  ex  coordinate  direction  this  auxiliary  system  takes 
the  form: 


d&  d<&  . 

Ao^t+Ay^  +  Az~dF  +  ^-°’ 


where  $  is  a  vector  of  length  rn(p  +  1)  for  a  boundary  condition  order  of  p  and  a  hyperbolic  system 
with  m  independent  variables.  Here  one  set  of  m  variables  is  the  trace  of  the  interior  solution. 
To  close  this  system  boundary  conditions  are  required  at  face  edges.  These  in  turn  are  provided 
by  edge  variables  which  satisfy  their  own  auxiliary  system.  Finally  the  edge  variables  require 
boundary  conditions  at  corners,  which  are  provided  by  another  set  of  auxiliary  corner  variables. 
Stein  successfully  implemented  all  of  these  coupled  systems.  His  code  achieves  the  a  priori  error 
estimates  of  [26].  We  illustrate  this  in  Figures  1-2.  In  the  first  case,  a  duct,  the  auxiliary  system  is 
closed  using  a  physical  boundary.  In  the  second,  free  space,  it  is  closed  using  the  edge  and  corner 
systems. 


2.2  Double  absorbing  boundary  implementation  of  CRBCs 

There  are  some  restrictions  and  inconveniences  with  the  implementation  of  CRBCs  described  above. 
In  particular  the  formulation  requires  writing  the  system  in  terms  of  characteristic  variables  as  well 
as  the  effective  use  of  sparse  matrix  solvers  (we  have  used  SuperLU  and  Umfpack)  to  treat  the 
corner  and  edge  systems.  The  former  issue  is  particularly  felt  when  solving  wave  equations  in 
second  order  form  or  when  using  the  popular  Yee  scheme  for  Maxwell’s  equations.  To  deal  with 
these  issues  we  have  developed  an  alternative  formulation  called  the  double  absorbing  boundary 
(DAB)  [23,  9].  In  the  DAB  we  solve  the  auxiliary  system  for  $  not  on  the  boundary  but  in  a  thin 
(one  element  or  one  stencil  width)  layer  near  the  boundary.  The  DAB  does  not  require  first  order 
systems,  characteristic  variables,  or  semi-implicit  systems  at  edges  and  corners.  It  is  slightly  more 
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3d  Duct  Waveguide  using  Complete  Radiation  Boundary  Conditions 


Figure  1:  Maximum  observed  and  predicted  errors  for  solution  of  the  acoustic  system  with  c  =  1  in 
a  duct  with  square  cross-section.  Here  the  solution  was  produced  by  a  point  source  at  the  origin, 
and  simulated  using  8th  order  differencing  in  space  combined  with  4th  order  Runge-Kutta  time 
stepping.  The  duct  width  is  2  as  is  the  width  of  the  computational  domain. 


expensive  than  the  original  formulation  in  terms  of  flops,  but  requires  less  memory  when  edges  and 
corners  are  present.  It  is  being  used  in  part  of  the  software  implementation  of  CRBCs  discussed 
below. 

2.3  Extensions  to  elastic  waves  and  stratified  media 

The  theoretical  developments  and  implementations  outlined  above  assume,  in  the  far  field,  that  the 
medium  is  homogeneous  and  isotropic  and  that  there  is  a  single  wave  speed.  This  is  reasonable, 
for  example,  for  studying  electromagnetic  scattering  in  a  vacuum  or  a  dielectric  medium,  but  is 
not  accurate  for  many  problems  arising,  for  example,  in  seismic  wave  studies.  As  such  we  have 
begun  the  process  of  extending  the  CRBCs  to  these  more  complex  situations.  In  particular  we  have 
demonstrated  excellent  accuracy,  even  without  full  parameter  optimization,  for  applications  of  the 
CRBCs  to  problems  in  stratified  media  as  well  as  for  some  anisotropic  problems  [21,  7,  22].  This 
allows  their  application,  for  example,  in  layered  waveguides  and  for  advective  acoustics. 

We  have  also  initiated  their  study  for  applications  to  elastic  waves  [8,  35,  34],  Complications 
with  the  elastic  wave  system  include  the  presence  of  multiple  wave  speeds,  boundary  conditions 
involving  oblique  derivatives,  and  anisotropic  media.  In  many  cases  stability  is  difficult  to  obtain; 
for  example  no  known  stable  PMLs  exist  for  certain  orthotropic  media  [6].  We  have  shown  that 
CRBCs  applied  to  the  elastic  wave  equation,  with  appropriate  parameter  choices,  are  always  long- 


4 


Free  Radiating  Box,  Absolute  Errors  and  Reflection  Coefficients,  8x=2/502, 8t=1/1 250, 0  <  t  <  100,  x  e  [-1,1] 
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Figure  2:  Maximum  observed  and  predicted  errors  for  solution  of  the  acoustic  system  with  c  =  1 
for  a  solution  in  free  space  produced  by  a  point  source  at  the  origin,  and  simulated  using  8th 
order  differencing  in  space  combined  with  4th  order  Runge-Kutta  time  stepping.  The  width  of  the 
computational  domain  is  2. 


time  stable.  This  is  in  contrast  with  all  other  high-order  radiation  conditions  which  have  been 
proposed  for  elastic  waves  which  suffer  from  long-time  error  growth.  We  illustrate  this  in  Figure  3 
below. 

We  are  actively  working  on  improving  these  extansions.  The  fundamental  issue,  being  studied 
by  a  Ph.D.  student,  John  Lagrone,  who  has  been  funded  as  an  RA  on  this  project  as  well  as  on 
a  related  STTR,  is  the  computation  of  optimal  parameters  and  an  estimate  of  convergence  with 
increasing  boundary  condition  order.  This  is  important  not  only  for  elastic  wave  models  but  also 
for  electromagnetic  waves  in  more  complex  media  such  as  nretanraterials.  He  is  also  looking  into 
more  general  formulations  in  case  they  are  needed  to  produce  efficient  methods  for  arbitrary  wave 
systems. 

2.4  A  CRBC  Software  Library 

In  work  with  HyPer Comp,  primarily  funded  by  an  STTR  but  also  partially  supported  by  this  grant, 
we  have  developed  an  open  source  software  library,  rbcpack  (www.rbcpack.org),  with  implemen¬ 
tations  of  CRBCs  for  various  volume  discretizations  of  Maxwell’s  equations.  The  library  generally 
involves  a  preprocessor  which  computes  optimal  parameters  given  e,  c,  T,  and  5,  sets  up  the  auxil¬ 
iary  system,  and  provides  an  interface  to  the  volume  solver.  Currently  two  pieces  to  the  library  are 
complete  or  nearly  complete.  The  first  is  the  implementation  in  the  popular  Yee  scheme,  developed 
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Figure  3:  The  magnitude  of  a  solution  component  for  a  simulation  of  transversely  1-periodic  elastic 
waves  in  2  +  1  dimensions.  The  solution  displaying  long-time  instability  was  obtained  using  the 
traditional  implementation  of  high-order  radiation  boundary  conditions  as  first  proposed  by  Hig¬ 
don.  The  long-time  stable  solutions  were  produced  using  our  new  CRBC  formulation  but  different 
parametrizations . 


by  John  Lagrone  and  utilizing  the  DAB  formulation.  This  has  been  completely  tested  with  our 
own  codes  and  extensively  documented;  it  will  be  released  once  our  first  external  users  verify  that 
it  can  be  linked  easily  as  a  library.  The  second  is  a  boundary-based  implementation  for  arbitrary- 
order  discontinuous  Galerkin  discretizations.  This  has  been  completely  tested  for  hexahedral  grids 
and  for  tetrahedral  grids  in  waveguides.  Once  the  implementation  on  tetrahedral  grids  for  exterior 
problems  is  complete,  which  we  expect  soon,  this  code  will  also  be  released.  It  can  be  used  in 
conjunction  with  HyPerComp’s  HDPhysics  platform.  The  third  component  is  a  DAB  implementa¬ 
tion  for  high-order  difference  approximations  to  second  order  formulations  with  Overture’s  cgrnx 
code  as  a  target  application.  Currently  a  two-dimensional  version  of  the  code,  developed  by  Fritz 
Juhnke,  a  doctoral  student  who  was  supported  as  an  RA  on  this  project  and  the  related  STTR,  is 
being  coupled  with  cgmx;  he  has  successfully  tested  it  using  a  home-grown  version  of  the  6th  order 
space-time  approximation  native  to  cgmx.  Once  this  coupling  is  proven  to  work  he  will  move  on  to 
the  three-dimensional  implementation.  Note  that  with  a  DAB  formulation  extensions  from  two  to 
three  space  dimensions  do  not  require  conceptually  distinct  discrete  systems,  so  it  is  reasonable  to 
expect  that  the  extension  will  be  complete  in  less  than  a  year. 
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2.5  Time-domain  integral  equations  for  electromagnetic  scattering 

An  alternative  to  the  volume-based  methods  we  work  on  for  simulating  electromagnetic  waves  are 
methods  based  on  integral  equations.  These  have  the  advantage  that  no  volume  grid  is  required,  so 
that  scattering  from  extremely  complex  structures  (e.g.  circuit  boards)  can  be  carried  out  without 
the  need  for  a  costly  grid  generation  step.  This  approach  has  been  reinvigorated  in  the  past  decade 
by  the  introduction  of  the  so-called  Plane  Wave  Fast  Time  Domain  Algorithm  (PWFTD)  [16,  36], 
which  is  a  time-domain  version  of  Greengard  and  Rokhlin’s  celebrated  fast  multipole  algorithm. 
However,  the  theory  behind  these  equations  and  the  algorithms  used  to  solve  them  is  not  well 
developed,  and  a  number  of  seemingly  attractive  discretization  techniques  experience  unexplained 
instabilities.  Our  first  goal  is  to  develop  more  efficient  time-stepping  procedures  for  these  unusual 
equations.  By  studying  the  simple  problem  of  scattering  from  a  sphere  we  have  established  a 
connection  between  these  time-domain  equations  and  neutral  delay  differential  equations,  and  have 
developed  stable  explicit  methods.  In  addition  we  have  shown  how  to  alter  the  integral  equation 
so  that  the  decay  properties  of  the  potentials  match  those  of  the  physical  fields.  This  is  joint  work 
with  Leslie  Greengard  of  the  Simons  Foundation  and  the  Courant  Institute  and  Charlie  Epstein  of 
the  University  of  Pennsylvania  [15].  In  Figure  4  we  illustrate  the  remarkably  different  long  time 
behavior  of  the  density  for  different  formulations  of  the  same  scattering  problem.  Note  that  the 
solution  of  the  wave  equation  itself  decays  exponentially. 


Solution  with  nonoscillatory  data:  ct=p=0  Solution  with  nonoscillatory  data:  a=1 ,  p=1/2 


Figure  4:  Density  fi{t)  at  mode  0  for  scattering  of  a  plane  Gaussian  pulse  from  a  sphere.  On  the 
left  the  solution  using  the  standard  single  layer  potential  and  on  the  right  an  appropriately  chosen 
combined  field  representation.  We  use  an  explicit  multistep  method  of  order  6. 


In  addition  we  have  developed  codes  for  the  accurate  evaluation  of  exact  solutions  for  both  scalar 
wave  and  electromagnetic  scattering  from  a  sphere  [18,  19].  These  provide  benchmark  solutions  for 
any  time-domain  scattering  solver. 

2.6  Frequency  domain  applications 

Although  the  primary  focus  of  our  boundary  condition  development  is  the  time-domain,  the  CRBCs 
can  also  be  used  for  problems  posed  in  the  frequency  domain.  This  development  was  led  by  Dr. 
Seungil  Kim,  who  for  two  years  was  a  postdoctoral  fellow  supported  by  the  grant.  The  basic 
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mathematical  formulation  and  the  theory  of  finite  element  approximations  is  worked  out  in  [30,  31]. 
In  particular  in  [31]  we  are  able  to  provide  a  complete  mathematical  justification  for  the  edge  and 
corner  closure  conditions,  which  as  yet  are  unavailable  in  the  time  domain. 

Besides  the  domain  truncation  itself,  CRBCs  can  be  used  in  place  of  PMLs  to  develop  sweeping 
preconditioners,  as  first  proposed  by  Engquist  and  Ying  [14].  These  address  one  of  the  most 
important  issues  for  frequency  domain  compuations,  namely  the  efficient  solution  of  the  resulting 
algebraic  system.  In  particular  it  is  argued  in  [14]  and  elsewhere  that  sweeping  preconditioners  can 
deliver  optimal  or  near-optimal  convergence  rates.  Kim  has  shown  that  CRBCs  are  an  effective 
replacement  for  PML  in  this  context  [32],  We  note  that  the  CRBCs  have  also  been  similarly  used 
to  develop  improved  parabolized  models  of  jet  noise  by  Towne  and  Colonius  [39,  40]. 

The  criteria  for  optimizing  the  CRBC  parameters  for  use  in  a  sweeping  preconditioner  differ 
somewhat  from  those  in  the  radiation  boundary  condition  problem,  and  in  any  case  most  appli¬ 
cations  involve  variable  coefficients  where  optimal  parameters  are  unknown.  Lagrone  has  done 
some  work  to  use  rational  Krylov  methods  (e.g.  [33])  to  construct  locally  optimized  transmission 
conditions.  We  expect  to  return  to  this  development  in  the  future. 

3  Energy-stable  high-resolution  volume  discretizations 

Our  belief  is  that  to  best  exploit  advances  in  computing  technology,  algorithms  based  on  high-order 
space-time  discretizations  are  needed.  The  reasons  for  this  are  two-fold: 

i.  One  aims  to  solve  more  difficult  problems  as  measured  by  the  number  of  wavelengths  over  which 

accuracy  must  be  maintained,  and  in  the  treatment  of  heterogeneous  media  with  disparate 
physical  properties.  More  simply  put,  one  wishes  to  solve  problems  exhibiting  multiple  spa- 
tiotemporal  scales.  Due  to  their  minimal  dispersion  error,  high-order  methods  require  orders 
of  magnitude  fewer  degrees-of-freedom  than  traditional  lower  order  schemes  when  the  number 
of  wavelengths  propagated  becomes  large,  and  thus  they  are  a  natural  choice  for  challenging 
problems. 

ii.  Higher  order  methods  generally  involve  more  localized  computation  than  their  lower  order  coun¬ 

terparts.  Thus  they  exhibit  larger  computation-to-communication  ratios  and  can  therefore 
better  exploit  modern  multicore  architectures. 

The  main  challenges  to  realizing  the  potential  of  high-order  discretizations  are  to  develop  methods 
which  are  robustly  stable  and  capable  of  dealing  with  complex  geometric  features.  For  unstructured 
meshes,  upwinded  discontinuous  Galerkin  (DG)  methods  [27]  are  an  excellent  choice  which  we 
employ.  In  particular  they  are  energy-stable  and  applicable  on  very  general  meshes.  The  downside  of 
DG  methods,  or  more  generally  any  element-based  methods,  is  artifically  stiff  derivative  operators. 
(See,  e.g.,  [42].)  The  effect  of  this  artificial  stiffness  is  the  need  to  take  small  times  steps,  which  has 
essentially  limited  most  methods  to  less-than-optimal  orders.  For  polynomial-based  discretizations 
the  only  way  to  avoid  this  artificial  stiffness  is  to  avoid  differentiation  throughout  the  polynomial’s 
domain  of  definition,  which  is  best  achieved  by  using  structured  grids.  Thus  we  are  focused  on  the 
use  of  hybrid  structured-unstructured  grids,  combining  upwind  DG  discretizations  with  energy- 
stable  structured  grid  methods  which  allow  large  time  steps  throughout  much  of  the  domain. 
We  are  studying  two  distinct  methods  of  this  type:  novel  spectral  elements  based  on  Hermite 


interpolation  [17,  5],  and  compact  difference  methods  derived  from  DG  ideas  [10].  In  addition  we 
have  developed  new  upwind  DG  methods  for  wave  equations  in  second  order  form. 

3.1  Hermite  Methods 

Hermite  methods  are  spectral  element  methods  with  unique  properties  which  make  them  ideal  for 
efficient  implementation  on  modern  multicore  architectures,  a  project  underway  in  collaboration 
with  Warburton  and  his  group  [41].  In  comparison  with  existing  techniques,  they  allow  the  in¬ 
dependent  evolution  of  thousands  of  degrees-of-freedom  over  relatively  large  time  steps, 
effectively  minimizing  communication  requirements.  In  particular,  the  only  time  step  restriction 
is  the  physical  CFL  condition  independent  of  order.  Thus  for  methods  of  order  10  or  more  the 
method  is  at  least  10  times  faster  than  standard  DG,  and  has  much  better  (essentially  optimal) 
computation-to-communication  characteristics. 

Hermite  methods  are  energy-stable,  but  the  analysis  requires  the  introduction  of  an  unusual 
energy  based  on  high-order  derivatives.  Precisely,  for  methods  of  order  2m  +  1,  we  prove  stability 
using  the  seminorm: 


This  leads  to  theoretical  barriers  to  proving  the  stability  for  coupling  with  standard  DG  schemes 
on  hybrid  grids.  Nonetheless  we  do  have  experimental  verifications  of  the  stability  of  hybridized 
Hermite-DG  methods  [11].  As  an  example,  we  have  solved  the  TM  Maxwell  system  on  a  disk 
using  the  hybrid  grid  shown  in  Figure  5,  with  convergence  results  for  our  9th  order  implementation 
tabulated  in  Table  1.  Of  note,  besides  the  fact  that  convergence  at  design  order  is  observed,  is  the 
ratio  of  the  time  steps  -  in  the  extreme  case  we  take  80  steps  in  each  DG  cell  for  every  step  in  a 
Hermite  cell.  Analyzing  the  stability  of  the  hybrid  method  is  a  priority  for  future  research. 


Figure  5:  Hybrid  grids  used  for  a  coupled  Hermite-DG  solver  for  Maxwell’s  equations.  We  use  the 
DG  scheme  on  the  mapped  triangular  elements  near  the  PEC  boundary  and  the  Hermite  scheme 
elsewhere. 
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h 

K 

CFL 

Error 

rate 

0.1 

447 

0.7 

1.83(E-03) 

1.37(E-08) 

0.08 

666 

0.7 

8.56(E-04) 

1.42(E-09) 

10.1 

0.05 

1525 

0.7 

4.09(E-04) 

1.68(E-11) 

9.4 

Table  1:  Experimental  convergence  for  the  hybrid  Hermite-DG  Maxwell  solver  on  a  disk. 


In  a  series  of  papers  we  have  examined  the  basic  properties  of  Hermite  methods,  including 
their  performance  for  nonsmooth  solutions  [2]  and  their  dispersion/dissipation  characteristics  [28, 
29].  The  latter  work  was  a  central  component  of  Chang- Young  Jang’s  doctoral  thesis;  Dr.  Jang 
received  some  RA  support  from  the  grant  while  he  was  a  student.  We  note  that  the  performance 
of  the  methods  for  propagating  discontinuous  solutions  is  particularly  impressive.  Interestingly  it 
improves  with  increasing  order;  see  Figure  6  below  for  comparisons  with  a  WENO  scheme  [37].  We 
currently  have  no  theory  which  fully  explains  this  excellent  behavior,  but  hope  to  develop  one  in 
the  future.  We  also  note  that  versions  with  full  shock-capturing  capabilities  are  being  developed 
by  our  collaborator  Appclo  and  his  student  Kornelus. 

Additional  work  includes  the  implementation  of  a  p-refinement  strategy  [12]  and  the  outlines 
of  a  scheme  for  /i-refinement  on  octrees  [1].  From  the  theoretical  perspective  we  need  to  prove 
stability  for  the  p-refinements.  We  also  need  to  develop  an  effective  strategy  for  choosing  between 
the  two  approaches. 

Lastly,  in  collaboration  with  Daniel  Appclo  of  UNM,  we  are  about  to  release  a  basic  open-source 
library,  CHIDES  (www.chides.org),  to  help  explain  Hermite  methods  and  enable  their  wider  use. 

3.2  Galerkin  Difference  Methods 

High-order  difference  methods  provide  a  simple  alternative  to  Hermite  methods  as  a  structured-grid 
solver.  Recently,  methods  with  the  summation- by-parts  (SBP)  property  have  been  incorporated 
into  many  unsteady  aerodynamics  simulators  due  to  their  guaranteed  stability  on  poor  grids  [38]. 
However,  a  defect  of  existing  high-order  SBP  methods  is  severe  loss  of  accuracy  at  boundaries. 
Galerkin  difference  methods  simply  use  the  local  piecewise  polynomial  bases  associated  with  stan¬ 
dard  difference  schemes  in  a  Galerkin  framework.  They  are  closed  at  boundaries  either  using 
extrapolation  to  ghost  points  or  by  leaving  the  basis  functions  associated  with  ghost  points  in  the 
Galerkin  basis.  In  [10]  we  show  that  the  resulting  compact  difference  schemes  are  super- 
convergent,  like  their  DG  counterparts,  but  have  essentially  order-independent  time 
step  stability  restrictions  just  as  central  difference  methods  possess  when  there  are  no  bound¬ 
aries.  See,  for  example,  Table  2  for  the  extremely  slow  growth  with  p  of  the  eigenvalues  of  the 
differentiation  matrices  including  the  boundary  closures.  Additionally,  although  the  mass  matrices 
are  not  diagonal,  they  are  tensor  products  of  banded  matrices  and  thus  can  be  inverted  with  linear 
cost. 

As  they  are  based  on  the  Galerkin  framework,  stable  hybridization  with  standard  DG  methods 
on  hybrid  grids  is  automatic. 

The  Galerkin  difference  methods  have  both  advantages  and  disadvantages  relative  to  Hermite 
schemes.  On  standard  architectures  they  are  somewhat  more  efficient  in  terms  of  flops,  but  their 
communication  patterns  are  not  as  ideal.  We  are  just  starting  to  experiment  with  the  schemes,  so 
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V 

3 

5 

7 

9 

11 

13 

15 

17 

Ax  •  p{Dp'u) 

2.02 

2.17 

3.41 

3.18 

3.67 

3.81 

3.57 

3.84 

A x-  p{Dp’x) 

2.02 

2.17 

2.27 

2.33 

2.39 

2.43 

2.46 

2.49 

Table  2:  Table  of  maximum  eigenvalues  and  normalized  stability  constraint  for  first  derivatives 
with  the  ghost,  G,  and  extrapolation,  X,  basis  closures  at  various  orders.  Here  p  is  the  method 
design  order  and  Dp  is  the  first  derivative  matrix  including  boundary  conditions. 


the  development  of  more  serious  multidimensional  implementations  is  a  subject  for  future  work, 
being  carried  out  in  collaboration  with  Jeff  Banks  of  RPI.  In  particular  we  are  also  considering 
embedded  boundary  formulations  which  would  enable  the  solution  of  problems  in  complex  geometry 
on  purely  Cartesian  grids. 


3.3  Upwind  DG  methods  for  second  order  wave  equations 

In  many  cases  second  order  formulations  of  wave  equations  have  advantages  over  first  order  for¬ 
mulations.  However,  in  our  view  the  usual  DG  schemes  for  wave  equations  in  second  order  form 
such  as  SIPDG  [20]  and  LDG  [13]  are  less  attractive  than  their  first  order  counterparts  due  to  the 
lack  of  a  natural  upwinding  strategy.  In  joint  work  with  Appelo  [4]  we  propose  a  new,  and 
quite  general,  formulation  of  DG  methods  for  wave  equations  in  second  order  form, 
with  applications  to  elastic  waves  demonstrated  in  [3].  The  essential  ideas  behind  the  method  are: 

i.  Introduce  a  weak  approximation,  v,  to  the  time  derivative  of  the  solution,  u. 


ii.  Build  fluxes  based  on  the  Lagrangian  form. 


Precisely,  consider  a  system  with  a  potential  energy  function  G( u,  Vu,  x).  Then  the  weak  form  on 
an  element  Clj  is  given  by: 
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where  v*  and  w*k  are  appropriately  chosen  fluxes.  Although  the  formulation  may  look  complex, 
its  implementation  is  fairly  straightforward,  and  for  complex  systems  it  uses  significantly  fewer 
unknowns  than  LDG.  In  [4]  we  observe  optimal  convergence  for  both  dissipative  upwind  fluxes  and 
energy-conserving  alternating  fluxes.  We  also  prove  optimal  convergence  in  the  energy  norm  for 
simple  cases. 
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The  method  is  naturally  formulated  for  any  system  in  Lagrangian  form,  and  we  have  exper¬ 
imented  with  it  for  various  nonlinear  model  problems.  These  often  develop  singularities,  and  an 
interesting  open  problem  is  to  develop  conditions  on  discretization  schemes  which  guarantee  con¬ 
vergence  to  physically  relevant  weak  solutions.  We  also  note  that  the  new  formulation  can  be  used 
to  derive  Galerkin  difference  methods  which  can  be  applied  on  hybrid  grids.  These  issues  will  be  a 
focus  of  our  continued  research  into  these  methods. 

4  Applications 

Although  our  overall  concentration  is  on  the  development  of  general-purpose  algorithms  for  waves, 
we  are  interested  in  applications  to  challenging  problems.  One  example  is  the  simulation  of  seismic 
waves.  This  is  particularly  relevant  for  modeling  earthquakes,  but  is  also  important  for  problems 
related  to  imaging  underground  geology.  For  the  former  problem  we  are  collaborating  with  Jacobo 
Bielak  of  Carnegie  Mellon  University  and  Dan  Givoli  of  the  Technion.  We  have  both  Hermite  and 
DG  elastic  wave  codes,  and  a  Galerkin  difference  elastic  code  is  planned,  but  we  need  to  look  at 
their  hybridization  near  geologic  features  in  the  earth’s  interior.  Developing  this  code  and  applying 
it  to  basic  benchmark  problems  will  be  a  next  step.  We  note  that  similar  ideas  can  be  used  in  the 
context  of  nondestructive  evaluation. 

In  collaboration  with  Daniel  Appelo  of  the  University  of  New  Mexico  and  Tim  Colonius  of 
Caltech  we  have  developed  a  Hermite-based  compressible  Navier-Stokes  solver.  A  major  application, 
primarily  funded  by  the  NSF  but  benefiting  from  the  work  done  under  this  grant,  is  the  direct 
numerical  simulation  of  jet  noise  at  Reynolds  numbers  an  order  of  magnitude  higher  than  those 
attained  previously,  to  within  one  order  of  magnitude  of  the  Reynolds  numbers  where  most  relevant 
experiments  are  performed.  The  code  is  fully  parallelized  and  we  have  carried  out  a  number  of 
experiments  with  turbulent  compressible  flows  to  verify  its  accuracy.  As  we  complete  some  of  our 
higher  Reynolds  number  runs,  we  face  the  problem  of  extracting  useful  insights  from  the  results. 
Here  we  hope  to  leverage  work  in  model  reduction,  control  theory,  and  data  mining  to  use  our 
simulations  to  develop  better  models  of  noise  production. 

Lastly  we  note  that  we  expect  HyPerComp  to  make  extensive  use  of  our  radiation  boundary 
condition  library  to  solve  problems  in  electromagnetism  and  acoustics  of  direct  interest  to  the  DOD. 
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